#include "Data.h"

/*Funcao: NewMatrixA
  Paraetros:
  float **matrixA - matrix original
  float *lineA - linha usada para o teste apos o metodo
  Saida:
  VOID
  A funcao salva a linha de teste da matrixA e divide os elementos pelo seu elemento da diagonal correspondente
*/
void newMatrixA(float **matrixA, float *lineA){
  int line, col;
  /* Loops irao percorrer toda a matrix, dividindo cada elemento pelo respectivo elemento da diagonal
     O proprio elemento da diagonal nao eh modificado
     A linha para o teste final eh copiada para o vetor lineA
  */
  for(line=0; line<J_ORDER; line++){
    for(col=0; col<J_ORDER; col++){
      if(line == J_ROW_TEST){ // Se a linha atual eh a de teste, copia elemento para o vetor lineA
        lineA[col] = matrixA[line][col];
      }
      if(line != col){ // Se nao for um elemento da diagonal
        matrixA[line][col] /= matrixA[line][line];
      }
    }
  }
}

/*Funcao: newVectorB
  Parametros:
  float **matrixA
  float *vectorB
  float *vectorX
  float *rB
  Saida:
  VOID
  A funcao gera o novo vetorB, dividido pelo elemento de A, no qual é zerado ao fim da funcao. A copia do vetorB ao vetorX ocorre para dar inicio ao metodo 
*/
void newVectorB(float **matrixA, float *vectorB, float *vectorX, float *rB){
  int line;
  /*Divide o elemento do vetor B pelo seu respectivo elemento da diagonal na matrixA. Em seguida, o elemento da diagonal da matrix eh zerado e o vetor B eh copiado para a variavel float rB*/
  for(line=0; line<J_ORDER; line++){
    if(line == J_ROW_TEST){ // Se eh o elemento a ser testado ao fim
      (*rB) = vectorB[line];
    }
    vectorB[line] /= matrixA[line][line]; // atualiza valor de B[i]
    vectorX[line] = vectorB[line]; // copia para vetorX
    matrixA[line][line] = 0.0; //zera matrixA[i,i]
  }
}

/*Funcao: updateX
  Parametros:
  float *oldVectorX
  float *vectorX
  Saida:
  VOID
  Atualiza os vetores de X, de modo que o vetorX da iteracao atual comece zerado
*/
void updateX(float *oldVectorX, float *vectorX){
  int i;
  /*Loop percorre ambos os vetores, X e oldX salvando os resultados da iteracao anterior e preparando o vetor X para a atual*/
  for(i=0; i<J_ORDER; i++){
    oldVectorX[i] = vectorX[i];
    vectorX[i] = 0.0;
  }
}

/*Funcao: calculateXi
  Parametros: 
  float *vectorX 
  float **matrixA
  float *vectorB
  float *oldVector
  Saida:
  VOID
  Calcula o vetorX da iteracao atual
*/
void calculateXi(float *vectorX, float **matrixA, float *vectorB, float *oldVector){
  int index, col;
  /*Faz o produto vetorial da MatrixA com o vetorX da ultima iteracao e a subtracao com o vetor B.*/
  for(index=0; index<J_ORDER; index++){
    for(col=0; col<J_ORDER; col++){
      //Produto vetorial da matrixA com o antigo vetor X
      vectorX[index] += matrixA[index][col]*oldVector[col];
    }
    //Subtracao do vetor X com o vetor B
    vectorX[index] *= -1;
    vectorX[index] = (vectorB[index] + vectorX[index]);
  }
}

/*FUncao: calculateDiff
  Parametros:
  float *vectorDiff
  float *vectorX
  float *oldVector
  Saida:
  VOID
  Calcula a diferenca absoluta entre os vetores X de cada iteracao
*/
void calculateDiff(float *vectorDiff, float *vectorX, float *oldVector){
  int index;
  /*Loop para percorrer todas as posicoes dos vetores e calcular a diferenca absoluta de cada um*/
  for(index=0; index<J_ORDER; index++){
    vectorDiff[index] = fabs(vectorX[index] - oldVector[index]);
  }
}

/*FUncao: getMaxValueFromX
  Parametros:
  float *vectorX
  Saida:
  float Max
  Retorna o maior elemento do vetorX
*/    
float getMaxValueFromVector(float *vector){
  int i;
  float max = -1.0;
  /*Loop que percorre todo o vetor, comparando com o maior ja encontrado.*/
  for(i=0; i<J_ORDER;i++){
    if(max < fabs(vector[i])){
      max = fabs(vector[i]);
    }
  }
  return max;
}
    
/*FUncao: getMr
  Parametros:
  void
  Saida:
  float Mr -
  Retorna o maior elemento do vetorDiff
*/
float getMr(float *vectorDiff, float *vectorX){
  float mr;
  mr = getMaxValueFromVector(vectorDiff)/getMaxValueFromVector(vectorX);
  return mr;
}

/*FUncao: getResult
  Parametros:
  float *vectorDiff
  float *vectorX
  Saida:
  float result
  Produto escalar do vetor LineA e o VetorX
*/
float getResult(float *lineA, float *vectorX){
  int i;
  float result = 0.0;
  //Loop para calcular o valor aproximado do sistema de equacoes lineares
  for(i=0; i<J_ORDER; i++){
    result += lineA[i]*vectorX[i];
  }
  return result;
}

/*Funcao: Jacobi
  Parametros:
  float **matrixA
  float *vectorB
  float *lineA
  float *vectorX
  float *oldVectorX
  float *vectorDiff
  float *rB
  Saida:
  VOID
  Metodo de Jacobi, resolve o sistema de equacoes lineares proposto no arquivo de entrada
*/
void Jacobi(float **matrixA, float *vectorB, float *lineA, float *vectorX, float *oldVector, float *vectorDiff, float *rB){
  //Variaveis de controle do loop
  int iterations = 0;
  float e = 1.0;
  
  //Calcula a nova matrixA*
  newMatrixA(matrixA, lineA);

  //Calcula o novo vetorB* e inicializa X
  newVectorB(matrixA, vectorB, vectorX, rB);
  
  //Loop principal do metodo
  while(iterations < J_ITE_MAX && e >= J_ERROR){
    //Prepara o vetor X para o comeco da iteracao.
    updateX(oldVector, vectorX);
    
    //Calcula o vetor X da iteracao atual
    calculateXi(vectorX, matrixA, vectorB, oldVector);

    //Atualiza o vetor de diferenca entre o vetor X e o OldVector
    calculateDiff(vectorDiff, vectorX, oldVector);
    
    //Calcula o erro de truncamento
    e = getMr(vectorDiff, vectorX);

    //Incrementa a variavel de controle do numero de iteracoes
    iterations++;
  }
    
  //Metodo finalizado. O vetor X contem os valores aproximados das incognitas do sistema. Eh impresso entao a comparacao do valor obtido com o esperado(salvo no comeco do programa)
  printf("\n--------------------\n");
  printf("Iterations: %d\n", iterations);
  printf("RowTest: %d ", J_ROW_TEST);
  printf("=> [%f] =? %f\n", getResult(lineA, vectorX), *rB);
  printf("--------------------\n");
}

